/*
  ivl_pow.c  -  implementation of exponential functions

  This file is part of the IVL Library (LIBIVL).
  
  LIBIVL - library for interval computations 
  
  Copyright (c) 2007 Mateusz Paprocki

  The IVL Library is free software; you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published
  by the Free Software Foundation; either version 2.1 of the License, or
  (at your option) any later version.

  The IVL Library is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with the IVL Library; see the file LICENCE.  If not, write to
  the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
  Boston, MA 02110-1301, USA.
*/

#include "interval.h"

int ivl_pow(ivl_ptr ivl, ivl_srcptr x, ivl_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	mpfr_t x_left, tmp;

	if (DOM_UNDEFINED(x) || DOM_UNDEFINED(y))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (ivl_is_equal_u(y, 0)) {
			if (ivl_is_equal_u(x, 0) || (ivl_point_inf_p(x) != 0)) 
				ivl->domain = IVL_UNDEFINED;
			else
				ivl_set_u(ivl, 1);
		} else {
			if (mpfr_sgn(IVL_RIGHT(x)) < 0)
				ivl->domain = IVL_UNDEFINED;
			else {
				mpfr_init2(x_left, ivl_get_prec(x));
				
				if (mpfr_sgn(IVL_LEFT(x)) < 0) {
					mpfr_set_ui(x_left, 0, GMP_RNDD);
					ivl->domain = IVL_PARTIAL;
				} else {
					mpfr_set(x_left, IVL_LEFT(x), GMP_RNDD);
					ivl->domain = x->domain | y->domain;
				}
				
				if (mpfr_cmp_ui(x_left, 1) < 0) {
					if (mpfr_cmp_ui(IVL_RIGHT(x), 1) <= 0) {
						if (mpfr_sgn(IVL_RIGHT(y)) > 0) {
							if (mpfr_sgn(IVL_LEFT(y)) >= 0) {
								inexact_left = mpfr_pow(IVL_LEFT(ivl), x_left, IVL_RIGHT(y), GMP_RNDD);
								inexact_right = mpfr_pow(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDU);
							} else {
								inexact_left = mpfr_pow(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDD);
								inexact_right = mpfr_pow(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDU);
							}
						} else {
							inexact_left = mpfr_pow(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDD);
							inexact_right = mpfr_pow(IVL_RIGHT(ivl), x_left, IVL_LEFT(y), GMP_RNDU);
						}
					} else {
						if (mpfr_sgn(IVL_RIGHT(y)) > 0) {
							if (mpfr_sgn(IVL_LEFT(y)) >= 0) {
								inexact_left = mpfr_pow(IVL_LEFT(ivl), x_left, IVL_RIGHT(y), GMP_RNDD);
								inexact_right = mpfr_pow(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDU);	
							} else {
								mpfr_init2(tmp, ivl_get_prec(ivl));
								
								inexact_left = mpfr_pow(IVL_LEFT(ivl), x_left, IVL_RIGHT(y), GMP_RNDD);
								inexact_left |= mpfr_pow(tmp, IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDD);
								
								if (mpfr_nan_p(IVL_LEFT(ivl)) || mpfr_nan_p(tmp))
									ivl->domain = IVL_UNDEFINED;
								else {																		
									mpfr_min(IVL_LEFT(ivl), IVL_LEFT(ivl), tmp, GMP_RNDD);  
								
									inexact_right = mpfr_pow(IVL_RIGHT(ivl), x_left, IVL_LEFT(y), GMP_RNDU);												
									inexact_right |= mpfr_pow(tmp, IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDU);
									
									if (mpfr_nan_p(IVL_LEFT(ivl)) || mpfr_nan_p(tmp))
										ivl->domain = IVL_UNDEFINED;
									else						
										mpfr_max(IVL_RIGHT(ivl), IVL_RIGHT(ivl), tmp, GMP_RNDD);
								}
								
								mpfr_clear(tmp);
							}
						} else {
							inexact_left = mpfr_pow(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDD);
							inexact_right = mpfr_pow(IVL_RIGHT(ivl), x_left, IVL_LEFT(y), GMP_RNDU);
						}
					}
				} else {
					if (mpfr_sgn(IVL_RIGHT(y)) > 0) {
						if (mpfr_sgn(IVL_LEFT(y)) >= 0) {
							inexact_left = mpfr_pow(IVL_LEFT(ivl), x_left, IVL_LEFT(y), GMP_RNDD);
							inexact_right = mpfr_pow(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDU);
						} else {
							inexact_left = mpfr_pow(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDD);
							inexact_right = mpfr_pow(IVL_RIGHT(ivl), IVL_RIGHT(x), IVL_RIGHT(y), GMP_RNDU);
						}
					} else {
						inexact_left = mpfr_pow(IVL_LEFT(ivl), IVL_RIGHT(x), IVL_LEFT(y), GMP_RNDD);
						inexact_right = mpfr_pow(IVL_RIGHT(ivl), x_left, IVL_RIGHT(y), GMP_RNDU);
					}
				}
				
				mpfr_clear(x_left);
			}
			
			IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
		}
	}
	
	return inexact;
}

int ivl_pow_u(ivl_ptr ivl, ivl_srcptr x, unsigned long y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (y == 0) {
			if (ivl_is_equal_u(x, 0) || (ivl_point_inf_p(x) != 0)) 
				ivl->domain = IVL_UNDEFINED;
			else
				ivl_set_u(ivl, 1);
		} else {
			if ((mpfr_sgn(IVL_LEFT(x)) < 0) && (mpfr_sgn(IVL_RIGHT(x)) < 0)) {
				ivl->domain = IVL_UNDEFINED;
			} else {
				ivl->domain = (mpfr_sgn(IVL_LEFT(x)) < 0) ? IVL_PARTIAL : x->domain; 
				
				inexact_left = mpfr_pow_ui(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
				inexact_right = mpfr_pow_ui(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);
				
				IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
			}
		}
	}
	
	return inexact;
}

int ivl_pow_s(ivl_ptr ivl, ivl_srcptr x, long y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (y == 0) {
			if (ivl_is_equal_u(x, 0) || (ivl_point_inf_p(x) != 0)) 
				ivl->domain = IVL_UNDEFINED;
			else
				ivl_set_u(ivl, 1);
		} else {
			if ((mpfr_sgn(IVL_LEFT(x)) < 0) && (mpfr_sgn(IVL_RIGHT(x)) < 0)) {
				ivl->domain = IVL_UNDEFINED;
			} else {
				ivl->domain = (mpfr_sgn(IVL_LEFT(x)) < 0) ? IVL_PARTIAL : x->domain;
				
				if (y >= 0) {
					inexact_left = mpfr_pow_si(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
					inexact_right = mpfr_pow_si(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);
				} else {
					inexact_left = mpfr_pow_si(IVL_LEFT(ivl), IVL_RIGHT(x), y, GMP_RNDD);
					inexact_right = mpfr_pow_si(IVL_RIGHT(ivl), IVL_LEFT(x), y, GMP_RNDU);
				}
				
				IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
			}
		}
	}
	
	return inexact;
}

int ivl_pow_z(ivl_ptr ivl, ivl_srcptr x, mpz_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (mpz_sgn(y) == 0) {
			if (ivl_is_equal_u(x, 0) || (ivl_point_inf_p(x) != 0)) 
				ivl->domain = IVL_UNDEFINED;
			else
				ivl_set_u(ivl, 1);
		} else {
			if ((mpfr_sgn(IVL_LEFT(x)) < 0) && (mpfr_sgn(IVL_RIGHT(x)) < 0)) {
				ivl->domain = IVL_UNDEFINED;
			} else {
				ivl->domain = (mpfr_sgn(IVL_LEFT(x)) < 0) ? IVL_PARTIAL : x->domain;
				
				if (mpz_sgn(y) >= 0) {
					inexact_left = mpfr_pow_z(IVL_LEFT(ivl), IVL_LEFT(x), y, GMP_RNDD);
					inexact_right = mpfr_pow_z(IVL_RIGHT(ivl), IVL_RIGHT(x), y, GMP_RNDU);
				} else {
					inexact_left = mpfr_pow_z(IVL_LEFT(ivl), IVL_RIGHT(x), y, GMP_RNDD);
					inexact_right = mpfr_pow_z(IVL_RIGHT(ivl), IVL_LEFT(x), y, GMP_RNDU);
				}
				
				IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
			}
		}
	}
	
	return inexact;
}

int ivl_u_pow(ivl_ptr ivl, unsigned long x, ivl_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(y))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (x == 0) {
			if (ivl_is_equal_u(y, 0) || (ivl_point_inf_p(y) != 0)) 
				ivl->domain = IVL_UNDEFINED;
			else
				ivl_set_u(ivl, 0);
		} else {
			if (x <= 1) {
				inexact_left = mpfr_ui_pow(IVL_LEFT(ivl), x, IVL_RIGHT(y), GMP_RNDD);
				inexact_right = mpfr_ui_pow(IVL_RIGHT(ivl), x, IVL_LEFT(y), GMP_RNDU);
			} else {
				inexact_left = mpfr_ui_pow(IVL_LEFT(ivl), x, IVL_LEFT(y), GMP_RNDD);
				inexact_right = mpfr_ui_pow(IVL_RIGHT(ivl), x, IVL_RIGHT(y), GMP_RNDU);
			}
			
			IVL_CHECK(ivl, y->domain, inexact, inexact_left, inexact_right);
		}
	}

	return inexact;
}

